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I  Dlstr-'.siit  isi 


We  find  that  potential  energy  rather  than  wave  height  is  a  better  ex-j  AvBUaLii  ccfTe^ 
perimental  and  analytic  criterion  for  determining  when  wave  breaking!  r.nz/or — 

will  occur.  A  simple  two-dimensional,  periodic  algorithm  is  developed' Dist  I  Spociai 
and  used  to  compare  breaking  onset  criteria  for  energy  input  &om  (1)1  Ck  A 
converging  sidewalls,  (2)  a  submerged  disturbance  and  (3)  wave  focusing.j  ' 

Wave-breaking  criteria  (potential  energy  or  the  more  classical  peak-to-‘  _ 

peak  wave  height)  are  a  function  of  the  rate  of  energy  input.  Large 
plunging  waves  occur  for  large  energy  input  rates  with  a  smooth  transi¬ 
tion  to  smaller  spilling  waves  for  lesser  energy  input  rates.  The  first  two 
kinds  of  energy  input  show  similar  trends  in  the  limit  as  the  energy  input 
rate  becomes  small.  The  third  case,  wave  focusing,  is  the  subject  of  an  / 

ongoing  investigation.  The  effects  of  wave  modulation  and  reflection  are 
also  discussed.  _ '' 


1  Introduction 


The  overall  morphology  of  the  surface  hydrodynanucs  of  ship  wakes  is  strongly 
affected  by  breaking  waves,  especially  at  the  bow  and  in  the  near  wake.  There  are 
essentially  two  types  of  breaking  waves — plunging  breakers  (with  a  large  degree  of 
overturning)  and  spilling  breakers  (with  white  water  only  near  the  crest).  Plunging 
breakers  are  an  important  factor  in  the  overturning  of  ships  in  rough  seas,  and  they 
often  form  continuously  at  the  bow,  producing  bubbles  and  foam  that  strongly 
affect  the  signature  of  a  ship  wake.  Spilling  breakers  are  more  common  in  the 
open  ocean  (due  to  wind)  and  in  breaking  wave  experiments;  they  also  occur  in 
the  near-ship  Kelvin  wave  pattern. 

A  recent  discussion  of  ship  wake  hydrodynamics  and  the  related  remote  sensing 
issues  is  given  by  Griffin  et  al.  (1989)  and  extensive  summaries  of  breaking  waves 
are  available  (Cokelet,  1978;  Kjeldsen  and  Myrhaug,  1978;  Griffin,  1984).  While 
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the  mechanisms  for  plunging  and  spilling  breaking  waves  are  often  thought  to  be 
quite  different,  we  will  show  that  their  formation  is  similar  with  their  qualitative 
features  depending  mainly  on  the  energy  input  rate. 

The  fundamental  experiments  for  determining  two-dimensional  wave-breaking 
criteria  (apart  from  that  induced  by  wind-wave  interaction)  fall  into  three  main 
categories:  (1)  the  focusing  of  essentially  two-dimensional  waves  in  the  lateral  di¬ 
rection  (Ramberg  and  Griffin,  1987;  Van  Dom  and  Pazan,  1975);  (2)  the  towing  of 
a  submerged  object  such  as  a  hydrofoil  (Duncan,  1981,  1983);  and  (3)  the  focusing 
of  variable  length  waves  from  a  modulated  wavemaker  or  wave  source  (Dommer- 
muth  et  al.,  1988;  Duncan  et  al.,  1987).  These  experimental  studies  propose  a 
wave-breaking  criterion  based  on  the  peak-to-peak  (crest-to-trough)  wave  heights. 

However,  the  validity  of  this  standard  criterion  has  been  questioned  (Melville 
and  Rapp,  1988),  in  part  because  peak-to-peak  wave  heights  vary  significantly 
during  breaking  and  often  decrease  just  before  breaking.  Any  criterion  will  be 
complicated  since  wave  breaking  is  not  amenable  to  analysis,  and  experimentally 
or  computationally  determined  criteria  will  be  a  fimction  of  many  parameters. 
Extensive  discussions  of  breaking  criteria  baised  on  wave  height  are  given  in  Ochi 
and  Tsai  (1983)  and  Huang  et  al.  (1986).  Breaking  criteria  based  upon  crest 
acceleration  are  discussed  by  Srokosz  (1986)  and  Longuet-Higgins  (1985).  An  ex¬ 
perimental  determination  of  the  onset  of  breaking  is  also  difficult  without  detailed 
velocity  measurements  at  the  crest  (Melville  and  Rapp,  1988;  Van  Dom  and  Pazan, 
1975),  which  are  usually  not  available  and  are  difficult  to  obtain. 

Computational  studies  of  breaking  waves  usually  apply  a  point  pressure  distur- 
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bance  (Longuet-Higgins  and  Cokelet,  1976)  or  obtain  breaking  conditions  simply 
from  having  energetic  initial  conditions  (Vinje  and  Brevig,  1981).  While  many  new 
algorithms  have  been  developed  for  breaking  waves,  no  study  has  systematically 
examined  the  parameters  that  cause  wave  breaking.  For  example,  wave  breaking 
caused  by  a  modulated  wavemaker  has  been  verified  by  computations  (Dommer- 
muth,  et  al.,  1988),  but  these  were  so  expensive  that  only  one  experimental  event 
was  verified.  In  addition,  previous  computations  tend  to  show  plunging  waves 
instead  of  the  more  commonly  observed  spilling  breakers. 

In  this  study,  we  computationally  examine  the  steepening  and  breaking  of  deep 
water  waves  generated  by  the  experimental  methods  cited  above.  We  consider  only 
spatially  periodic  computations,  so  an  ad  hoc  energy  input  term  is  deduced  for 
the  convergent  wave  channel.  Although  the  periodic  boundary  conditions  preclude 
studying  the  wavemaker  problem,  we  briefly  examine  the  eflTect  of  wave  modulation 
by  examining  a  larger  computational  region  (more  than  one  primary  wavelength) 
as  in  Dold  and  Peregrine  (1986).  The  effect  of  reflection  frmn  a  beach  can  be 
considered  by  putting  a  small  standing  wave  component  in  the  initial  conditions. 
Finally,  to  compare  to  the  second  type  of  experiment  we  use  an  array  of  submerged 
disturbances  (dipoles  in  this  case)  to  preserve  spatial  periodicity. 

There  are  difficulties  in  interpreting  the  differences  caused  by  the  computationad 
spatial  periodicity  as  compared  to  temporal  periodicity  (for  the  unmodulated  wave- 
maker)  of  the  experiments.  This  also  affects  comparisons  of  temporal  vs.  spatial 
growth  when  we  model  a  convergent  channel.  Also,  experiments  continue  after 
breaking  occurs,  while  the  time-marching  computations  must  stop  at  the  first  oc- 
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currence  of  breaking  unless  an  ad  hoc  condition  is  used  to  model  the  turbulence 
and  air  entrainment.  There  is  evidence  that  accurate  spectral  computations  break 
down  sooner  (Huh  and  Schultz,  1989),  indicating  the  possible  formation  of  a  sin¬ 
gularity  amd  the  failure  of  potential  theory  before  the  wave  crest  approaches  the 
forward  face. 

In  this  report,  experimental  data  for  the  wave  height  only  are  analyzed  since 
these  are  the  traditional  measurements  and  are  easily  obtained  compared  to  veloc¬ 
ity  measurements.  This  data  is  used  to  show  that  the  potential  energy  reduces  the 
experimental  scatter  in  breaking  criteria  and  to  show  it  acts  as  a  better  predictor 
(rather  than  just  indicator)  of  breaking. 

In  section  2,  we  pose  the  problem  for  periodic  waves,  including  the  modelling  of 
the  growth  in  energy  and  the  effect  of  submerged  disturbances.  Section  3  contains 
a  summary  of  earlier  computational  progress  followed  by  a  formulation  of  the  two 
computational  techniques  used  in  this  study.  Section  4  presents  numerical  results, 
including  comparisons  to  related  computationzd  schemes  and  the  development  of 
breaking  criteria  as  a  function  of  the  energy  input  rate.  The  computations  are 
compared  to  previous  experiments  in  section  5.  Measurements  of  momentum  and 
energy  losses  after  breaking  are  presented  in  section  6,  and  section  7  summarizes 
the  findings. 
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2  Problem  Formulation 


The  problem  domain  is  shown  in  Fig.  1.  The  scales  are  chosen  to  make  gravity  and 
the  prim2jy  wavenumber  unity.  The  phase  speed  and  angular  velocity  of  a  linear 
wave  will  then  be  imity  as  well.  The  initial  boundary  value  problem  solution  is 
described  by  a  complex  potential  to(^)  =  ^  +  where  ^  is  the  velocity  potential, 
is  the  stream  function  and  ^  =  2  +  ty  represents  the  two  spatial  coordinates. 
At  every  time  step,  the  unknown  boundary  values  of  the  velocity  potential  (half  of 
the  values  are  known  from  the  boundairy  conditions)  are  solved  using  the  Cauchy 
integral  theorem: 


where  a  is  0  or  2v  if  the  location  of  the  kernel  singularity,  Cib»  is  outside  or  inside 
the  boundary,  respectively.  If  the  kernel  singularity  is  on  the  boundary  (C*  €  Oil), 
a  is  equal  to  the  included  angle,  and  the  integral  is  treated  as  principal-valued. 

The  kinematic  and  dynamic  boimdary  conditions  of  a  free  surface  for  an  inviscid 
flow  are  given  as 


_  dw* 


(2) 


and 


D<f>  1  dw* 

~Di  ^  2 

Here,  p  is  a  prescribed  pressure  (normally  0,  as  in  this  study,  unless  surface  tension 
or  wind  effects  are  included),  DfDt  is  a  material  derivative,  and  *  denotes  the 
complex  conjugate.  The  kinematic  condition  requires  that  material  particles  on 
the  free  surface  stay  on  the  free  surface,  and  the  dynamic  boundao’y  condition 
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requires  that  the  pressure  remains  constant  at  the  free  surface  in  the  absence  of 
surface  tension. 

We  complete  the  problem  by  applying  a  periodic  boundary  condition  in  the 
horizontal  direction,  uj(^)  =  w{^  +  2x),  the  deep  water  condition  requiring  w  —*  0 
as  1/  —oo,  and  initial  conditions.  These  initial  conditions  on  the  free  surface 
can  be  homogene<~  us  if  the  waves  are  forced  by  submerged  dipoles.  For  purposes  of 
illustration  and  comparison,  we  generally  use  the  same  initial  conditions  as  Mclver 
and  Peregrine  (1981): 

Tf  =  a  sin  X  and  <^  =  a  cos  x.  (4) 

These  initial  conditions  satisfy  linear  theory  as  the  wave  amplitude  a  approaches 
0.  We  apply  increasing  values  of  a  until  the  wave  breaks.  We  also  apply  more  com¬ 
plicated  initial  conditions  to  examine  wave  modulation  and  the  effect  of  applying 
initial  conditions  computed  from  steadily  progressing  waveforms  using  the  method 
of  Schwartz  and  Vanden-Broeck  (1979). 

2.1  Modelling  the  convergent  wave  channel 

The  convergent  channel  is  inherently  a  three-dimensional  problem.  If  the  conver¬ 
gence  is  small  in  z,  the  spanwise  direction  out  of  the  plane  of  Fig.  1,  a  multiple- 
scales  approach  will  lead  to  the  three-dimensional  effect  being  delayed  to  a  Pois¬ 
son  equation  at  higher  order,  with  the  lowest-order  solution  being  that  of  a  non- 
convergent  wave  channel.  The  sidewall  boundary  condition  simply  becomes  = 
e<^x,  where  c  is  the  slope  of  the  converging  walls  and  the  subscripts  x  and  z  rep- 
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resent  partial  differentiation.  This  method  of  analysis  precludes  the  possibility  of 
modelling  the  flow  as  spatially  periodic  and  hence  makes  the  problem  intractable 
using  periodic  algorithms. 

If  instead  we  model  the  channel  as  one  with  straight  walls  converging  in  time, 
we  can  continue  to  model  the  problem  as  spatially  periodic.  We  assume  that  the 
convergent  channel  sidewall  boundary  conditions  can  be  approximated  by 


<f»z  =  ±ccj  at  z  =  z^  =  Zo±  ecgt. 


(5) 


where  Zo  is  the  scaled  width  of  the  channel  at  the  wavemaker  and  c-  is  the  av¬ 


erage  group  velocity  of  the  fundamental  wave.  These  conditions  lead  to  a  two- 
dimensional  Poisson  equation  with  a  small  right-hand  side: 


<l>xx  "b  4^yy  — 


tc„ 


Zo  —  ec„t 


(6) 


This  assumption  is  not  satisfactory  because  it  requires  significant  flow  at  infinite 
depths  and  consequently  the  Poisson  equation  cannot  be  made  into  the  easily 
computed  Laplace  equation. 

Instead  of  using  the  multiple-scales  approach  or  allowing  the  sidewalls  to  con¬ 
verge  in  time,  we  simply  add  a  term  in  the  Bernoulli  equation  that  causes  the 
energy  to  increase  exponentially  in  time  so  that  (3)  becomes 


D<i>  1 

Dt~  ^  2 


dw* 


+  l4>- 


(7) 


This  is  similar  (but  with  the  opposite  sign)  to  the  dissipative  effect  used  by  Lamb  to 
develop  the  radiation  condition.  Starting  with  a  wave  of  small  amplitude,  this  term 
eventually  causes  the  wave  to  break,  much  like  the  convergent  channel.  While  this 
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causes  exponential  temporal  growth  in  the  computations,  the  ex]}erimental  energy 
per  unit  width,  E,  grows  in  space.  The  corresponding  temporal  growth  parameter 
can  be  approximately  related  through  the  group  velocity  and  conservation  of  energy 


by 


1  dE 
E  dt  ~ 


1  dz„ 


dt  Zyj 


1  dz^  ec, 


a 

Xw' 


(8) 


The  dimensionless  values  (c  =  1/16)  for  the  experiments  considered  here  corre¬ 
spond  to 


7  =  a: 


II 

W' 


(9) 


where  the  dimensional  constant  K  is  .008m/s^,  and  W  and  T  are  the  dimensional 
values  of  the  local  channel  width  and  wavemaker  period,  respectively.  This  keeps 
the  growth  parameter  below  0.05  in  the  region  of  breaking. 


2.2  Modelling  periodic  submerged  disturbances 

We  wish  to  show  how  submerged  disturbances  can  force  wave  breaking  because  this 
case  is  far  easier  to  compute  accurately  than  surface-piercing  bodies.  To  keep  the 
periodic  boundary  condition  requires  that  the  disturbances  be  periodic.  Rather 
than  modelling  a  complex  two-dimensional  shape  such  as  a  hydrofoil,  we  use  a 
periodic  array  of  moving  dipoles; 

Wdp  =  arfpCot(i(^  -  ,  (10) 

where  a^p  and  (dp  are  the  strength  and  location  of  the  dipole,  respectively.  The 
dipole  depth,  ddp,  and  velocity,  Vdp,  are  prescribed  such  that  (dp  =  (vdpt,  —ddp).  As 
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long  as  the  strength  of  the  dipole  is  not  too  large  and  the  depth  is  not  too  small, 
the  dipole  can  closely  represent  a  cylinder  of  radius  r^p  if  =  Vipr\j2. 

Since  the  dipole  represents  a  simple  pole  in  the  complex  potential  plane,  the 
only  alteration  required  is  to  subtract  out  the  singular  pairt  and  solve  only  for  the 
remaining  regiilar  term  from  the  integral  equation  (1). 

3  Computations  of  Steep  and  Breaking  Waves 

3.1  Recent  computational  advances 

Although  formal  analytical  techniques  have  been  developed  for  small-amplitude 
gravity  waves,  unsteady  auid  steep  waves  must  be  solved  numerically.  The  most 
efficient  of  these  algorithms  are  based  on  boundary  integral  techniques.  Even  then, 
the  algorithms  can  be  rather  time  consuming.  Hence,  no  thorough  and  complete 
parametric  study  has  been  performed  on  gravity  waves.  Even  more  important,  to 
reduce  the  computational  effort,  the  problem  domadn  is  kept  as  small  as  possible 
by  applying  periodic  boundary  conditions.  Recently,  computations  with  many 
fundamental  wavelengths  inside  the  periodic  domain  have  been  applied  by  Dold 
and  Peregrine  (1986),  and  the  nonperiodic  fully  nonlinear  wavemaker  problem  has 
been  computed  by  Doiiunermuth  et  al.  (1938).  Casual  observations  of  breaking 
waves  show  that  they  are  not  spati2dly  periodic.  In  this  report,  we  only  present 
results  for  the  periodic  problem,  although  by  using  a  large  spatial  period,  the  model 
can  represent  results  on  an  infinite  domain. 

The  boundary  integral  numerical  schemes  for  irrotational  flow  problems  can  be 
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broken  into  three  general  approaches  based  on  Green’s  functions  (Longuet-Higgins 
and  Cokelet,  1976;  Vanden-Broeck,  1980),  vortex  dynamics  (Baker,  Meiron  and 
Orszag,  1982),  or  the  Cauchy  integral  theorem  for  complex  potentials  (Vinje  an'^' 
Brevig,  1981).  To  some  extent,  the  three  techniques  give  equivalent  results  (Mclver 
and  Peregrine,  1981).  Recent  work  (Bold  and  Peregrine,  1984)  has  shown  that  al¬ 
gorithms  based  on  the  Cauchy  Integral  theorem  can  be  up  to  50  times  faster  than 
Green’s  function  algorithms  and  10  times  faster  than  those  using  vortex  meth¬ 
ods.  Lin  et  al.  (1984)  use  the  Cauchy  formulation  when  solving  two-dimensional 
problems  and  revert  to  the  Green’s  function  algorithm  for  axisymmetric  problems. 
Apparently,  the  efficiency  of  the  complex  algebra  is  significant  amd  normal  deriva¬ 
tives  of  the  Green’s  function  need  not  be  found. 

Here  we  report  on  two  algorithms  based  on  the  Cauchy  integral  theorem.  The 
first  is  an  improvement  of  a  piecewise-linear  algorithm  of  Vinje  and  Brevig  (1981) 
as  described  in  Schultz  and  Hong  (1989).  The  second  is  a  spectral  technique  similar 
to  that  proposed  by  Roberts  (1983)  and  described  in  Huh  and  Schultz  (1989).  We 
describe  both  methods  because  the  first  method,  although  less  accurate  and  more 
computationally  intensive,  is  also  more  robust.  Comparisons  will  be  presented. 

In  both  methods,  the  physical  domain  is  mapped  to  an  approximate  unit  circle 
using  the  conformal  transformation: 

C  =  e"  (11) 

(see  Fig.  1).  This  eliminates  the  periodic  boundary  conditions  and  sharp  compu¬ 
tational  corners  used  by  Vinje  and  Brevig  (1981).  All  derivatives  are  taken  in  the 
conformed  space — the  piecewise-linear  method  uses  three-point  central  differences 
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while  the  spectral  method  tzdces  derivatives  in  the  spectral  space  of  the  conformed 
representation.  An  alternate  method  (not  used  here)  does  not  use  conformal  map¬ 
ping  but  replaces  the  infinitely  periodic  integrand  with  a  summation  over  a  finite 
domain  to  form  a  cotangent  kernel,  as  in  Baker  et  al.  (1982). 

The  algebraic  system  that  results  from  discretizing  the  integral  equation  is  it¬ 
eratively  solved  for  both  methods  using  a  generalized  minimum  residual  method 
(GMRES)  based  on  the  work  of  Saad  and  Schultz  (1986).  This  variation  of  the 
conjugate  gradient  method  for  nonsymmetric  matrices  works  very  well  on  the  diag¬ 
onally  dominant  matrices  of  either  method-especially  the  matrices  from  the  spec¬ 
tral  algorithm.  The  time  marching  is  also  similar  in  both  algorithms.  Fourth-order 
Runge-Kutta-Gill  and  Hamming  predictor-corrector  methods  with  an  automatic 
adjiistment  of  step  size  were  both  used,  with  the  predictor-corrector  method  show¬ 
ing  the  greater  computational  eflSciency,  especially  for  the  higher-accuracy  com¬ 
putations.  Filtering  of  the  spectral  computations  is  probably  necessary  to  allow 
them  to  proceed  closer  to  breaking,  although  none  vfas  performed  here.  Filtering 
is  discussed  in  Huh  and  Schultz  (1989). 

3.2  Piecewise-linear  computational  technique 

We  take  in  (1)  to  approach  the  boundary  from  the  outside  of  the  domain  so  that 
a  is  zero,  although  there  are  computational  reasons  for  placing  Ck  slightly  away 
from  the  contour  in  some  cases  (Schultz  and  Hong,  1989).  The  algebradc  system 
is  formed  by  discretization  of  (1),  as  explained  in  Vinje  and  Brevig  (1981),  and  by 
letting  the  kernel  singularity  approach  each  of  the  N  nodal  points,  i.e.,  (k  —*  (k-  A 
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speci2Ll  limiting  process  is  needed  to  evaluate  the  integration  near  Ck-  The  system 

of  linear  algebraic  equations  when  w  is  discretized  as  a  piecewise-linear  function 

between  the  N  boundary  nodes  is  as  follows: 

N 

=  0  for  =  (12) 

i=i 

where 

r,-,  =  InC^y-'  ~^)  -  )  for  j  #  k,  (13) 

Cj+1  -  —  Qk  —  €i-i  Ci-1  -  Qk 

Til  =  (14) 

Elq.  (13)  is  evaluated  using  I’Hopital’s  rule  when  j  i  +  1  or  —  1,  Moving 
the  known  boundary  conditions  to  the  right-hand  side  gives  a  complex  algebraic 
system  for  unknown  V’j  oii  the  free  surface. 

When  the  complex  potential  is  known  along  the  domain  boundary,  the  solution 
can  be  stepped  forward  in  time  using  the  Bemoulh  equation  and  the  kinematic 
boundary  conditions.  We  solve  this  problem  in  a  similar  way  to  Vinje  and  Brevig 
with  (1981)the  following  changes: 

1)  Rather  than  using  the  real  or  imaginary  peuts  of  the  discretized  Cauchy  inte¬ 
gral  theorem  (depending  on  whether  the  real  or  imaginary  part  of  w  is  known),  we 
use  both  to  give  2N  real  equations  and  N  real  unknowns.  Numerical  experiments 
for  known  test  cases  (Schultz  and  Hong,  1989)  show  that  the  least-squares  solution 
is  better  for  nearly  circular  contours,  especially  whai  the  node  placement  is  irreg¬ 
ular  (as  will  be  the  case  after  nodes  are  converted  on  the  free  surface);  however, 
both  results  are  second-order  convergent.  The  solution  time  for  a  direct  method  of 
inverting  the  overdetermined  system  would  be  twice  that  of  the  determined  system 
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but  our  experience  with  iterative  conjugate  gradient  solvers  indicates  an  increase 
in  computational  costs  of  only  10  percent. 

2)  We  use  a  conformal  map  to  eliminate  the  bottom  and  periodic  boundary 
conditions. 

3)  We  use  a  central  difference  form  for  dw/d^  (or  dw/d(),  while  Vinje  and  Bre- 
vig  use  a  truncated  analytic  form.  Since  the  solution  is  piecewise-linear  rather 
than  analytic,  we  have  found  that  some  numerical  instabilities  can  develop  using 
the  truncated  analytic  form.  One  can  easily  find  examples  where  the  derivative 
dwfd^  at  a  corner  of  the  contour  computed  using  the  analjrtic  form  lies  outside 
the  range  computed  by  the  forward  and  backward  derivatives.  This  violates  the 
spirit  of  using  piecewise-linear  functions  and  can  lead  to  numerical  instabilities, 
although  the  truncated  analjrtic  form  works  better  when  the  contour  is  smooth. 

3.3  Spectral  computational  technique 

Roberts  (1983)  used  a  desingularized  kernel  in  his  vortex  formulation.  Generally,  it 
is  difficult  to  find  a  suitable  desingularized  form  of  a  kernel  in  an  integral  equation, 
but  in  the  complex  formulation  it  is  relatively  simple.  The  Cauchy  integral  equation 
(3)  can  easily  be  rewritten  as 

where  the  principal  value  integral  can  be  replaced  by  the  closed  contour  integral 
since  the  integrand  is  no  longer  singular.  When  (  approaches  Cki  the  integrand 
approaches  dwfda  at  the  ibth  node  point.  Therefore,  this  kernel  does  not  show  the 
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singular  behavior  as  (  approaches  (k.  The  integral  equation  (1)  is  converted  to  the 


following  sets  of  equations  for  ib  =  1  ...  TV: 


s 

Elf  =  0 


for  A:  =  1, . . . ,  TV  , 


where  TV  is  the  number  of  nodes,  and  is  represented  by 


j.  ^1 


I  (ft.  ai  =  t. 

The  algebraic  system  (7)  effectively  becomes  a  differential  system  because  Ijk  in¬ 
cludes  the  derivative  of  to.  To  evaluate  these  derivatives  spectrally,  we  use  a  car¬ 
dinal  function  representation  of  to  (Boyd,  1989): 


where 


and  the  derivative  of  Cj  is 


Then,  (7)  becomes 


cot  f  (si  -  3^)  if  t  7^  j 
0  if  t  =  j  . 


J)rjjktOj  =  0  for  A:  =  1,...,TV  , 

i=i 


where  the  influence  coefficients  Fj*  are  now 


r,t  = 


.  “  Dill../*  7  =  * 
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Unlike  the  method  of  Baker  et  al.  (1982)  which  evaluates  the  integrand  at  every 
other  point,  the  desingularized  kernel  is  evaluated  at  every  nodal  point;  hence,  the 
matrix  is  full. 


4  Typical  Computational  Results 

4.1  Convergence  and  stability 

The  convergence  of  the  spectraJ  method  is  compared  to  the  piecewise-linear  method 
in  Huh  and  Schultz  (1989)  for  steady  test  problems.  Here,  we  briefly  extend  this 
comparison  to  time-marching  by  examining  a  gravity- wave  problem  with  the  initial 
conditions  given  by  (4)  and  a  =  0.2.  Fig.  2  compares  errors  using  IV  =  16  for  the 
conservation  of  energy — the  conservation  of  mass  results  are  very  similar.  The 
conservation  of  energy  is  determined  by  (E{t)  —  E{0))/E{0),  where  E  is  given  by 

(23) 

Not  only  are  the  resxilts  for  the  spectral  computations  much  more  accurate,  but 
they  are  less  computationally  intensive.  Ail  these  computations  for  a  =  0.2  can 
be  computed  indeflnitely  without  further  loss  in  accuracy  (as  tested  to  t  =  200). 
The  low  resolution  time  computations  were  made  with  Ct  =  Cc  =  10“^  (error 
criteria  for  the  time-mau’ching  and  iterative  algebraic  solver,  respectively),  the  high 
resolution  computations  with  c*  =  Cc  =  10~'°.  Further  reflning  the  time-marching 
parameters  for  either  method  does  not  improve  the  computational  error.  For  this 
initial  condition,  the  spectral  computations  can  give  essentially  double  precision 
(16  digit)  machine  accuracy  when  N  =  Z2  and  Cj  =  Cc  =  10“^^. 
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4.2  Steadily  progressing  wave 

We  first  show  that  under  special  circumstances  waves  of  large  amplitude  do  not 
break.  Specifically,  we  examine  gravity  waves  of  permanent  form  and  suppress  the 
Benjamin-Feir  instability  by  applying  periodic  boimdary  conditions  that  do  not 
allow  subharmonic  disturbances  (Longuet-Higgins,  1978).  The  initial  conditions  for 
a  steadily  progressing  wave  can  be  computed  from  a  series  expansion  as  performed 
by  Stokes  (1880)  and  extended  using  computer  Jilgebra  by  Schwartz  (1974).  Rather 
than  use  series  acceleration  techniques,  we  compute  the  initial  conditions  for  our 
time-marching  code  from  the  iterative  method  of  Schwartz  and  Vanden-Broeck 
(1979).  We  set  their  surface  tension  parameto’  to  zero  and  modify  their  mapping 
slightly  to  desingularize  the  mapping  at  very  high  amplitudes. 

Because  obtaining  accurate  initial  conditions  is  the  "weak  link”  in  these  compu¬ 
tations,  we  compute  twice  as  many  points  as  we  use  in  the  time-marching  algorithm 
and  discard  every  other  value.  To  obtain  an  accurate  Jacobian  matrix  for  proper 
convergence,  these  computations  must  be  performed  in  double  precision  (16  digits). 
Surprisingly,  unless  the  amplitude  is  very  near  the  limiting  Stokes  wave  height  that 
forms  the  120*’  crest,  single  precision  is  suffident  for  the  time-marching  program. 

Figs.  3a-b  show  eight  wave  height  profiles  spaced  At  =  0.5  apart  for  the  ampli¬ 
tude  parameter  (»7mo*— »7m»ti)/2x  =  0.1  and  0.115,  respectively.  The  time-marching 
computations  were  performed  using  N  =  64.  The  spectral  algorithm  was  used  in 
Fig.  3a,  while  the  more  robust  piecewise-linear  code  was  required  for  Fig.  3b  since 
the  amplitude  is  near  the  Stokes  limit  of  (t^ma*  ~  ’7min)/27r  »  0.14.  The  computed 
phase  velocities  (which  are  approximately  10.3  and  14.0  percent  higher  than  linear 
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waves  for  the  cases  in  Figs.  3a  and  3b,  respectively),  as  determined  by  marching  for 
long  time,  are  accurate  to  within  .005  percent.  The  RMS  variation  in  the  trough  to 
crest  wave  height  is  8x10“®  for  Fig.  3a  and  3x10“^  for  Fig.  3b.  The  overdetermined 
system  in  this  case  is  twice  as  accurate  as  the  strong  system  and  only  requires  20 
percent  more  CPU  time.  The  spectrad  computation  was  again  the  most  accurate, 
conserving  energy  and  mass  to  within  one  part  in  10®  (compatred  to  one  part  in 
10^  for  the  equivalent  piecewise- linear  computations  and  one  part  in  10^  for  the 
higher  wave  of  Fig.  3b).  The  potential  energy  was  constant  to  seven  significant 
digits  (PE  =  0.1440120,  KE  =  0.15158)  for  Fig.  3a  and  three  significant  digits  for 
Fig.  3b  (PE  =  0.1815,  KE  =  0.19435). 

These  waves  are  not  foimd  experimentally,  because  even  if  they  could  be  formed, 
they  are  subject  to  a  Benjamin-Feir  instability  (which  we  suppress  with  our  peri¬ 
odic  boundary  conditions). 

4.3  Simple  harmonic  linear  wave  initial  conditions 

Fig.  4  shows  typical  results  of  free  surface  profiles  for  dimensionless  times  0.1  apau't 
for  two  different  initial  amplitudes.  The  first  family  of  curves  for  a  =  0.3  (4)  results 
in  a  “spilling”  breaker.  When  more  energetic  initial  conditions  are  used,  as  in 
Fig.  4b  (a  =  0.544),  the  wave  becomes  a  “plimging”  breaker.  The  algorithm  breaks 
down  at  the  last  time  step  shown  because  of  insufficient  spatial  and/or  temporal 
resolution.  This  breakdown  exhibits  itself  as  a  fjulure  of  the  iterative  solution  of 
the  algebr2uc  system  to  converge.  The  wave  profiles  change  from  solid  to  dashed 
lines  when  the  spectral  solution  no  longer  converges.  Until  this  time,  the  spectral 


(solid  line)  and  piecewise-linear  (dashed  line)  computations  are  nearly  identical 
although  the  spectral  result  is  far  more  accurate.  (Slight  differences  in  the  profile 
are  shown  for  the  last  plotted  spectral  computation.)  This  is  an  indication  that  the 
spectral  computations  are  not  as  robust  as  the  piecewise-linear  computations.  Our 
experiences  have  shown,  however,  that  more  precise  piecewise-linear  computations 
(by  increasing  the  number  of  nodes  or  decreasing  the  time  step)  break  down  sooner. 
The  opposite  is  true  for  the  spectral  computations.  We  are  continuing  to  study  this 
in  an  attempt  to  see  if  singularities  develop  in  the  inviscid  model  before  “breaking” 
occurs,  as  in  the  study  of  singularity  formation  in  vortex  sheets  (Krasny,  1986). 

We  have  run  numerical  simulations  for  many  values  of  a  to  determine  the  initial 
conditions  (4)  that  cause  breaking  or  spilling.  We  find  that  for  a  >  0.28  waves  will 
spill  and  for  a  <  0.27  waves  will  progress  indefinitely.  Typical  computations  use 
iV  =  60  or  80  and  At  =  0.1  or  0.05.  These  results  are  somewhat  sensitive  to  the 
initisd  conditions  in  that  an  initial  condition  using  a  three-term  Stokes  profile, 

1/  =  a  sin  X  -I-  0.5a^  sin  2x  -f-  0.375a^  sin  3x,  (24) 

does  not  apply  as  large  a  perturbation  to  the  steady  form  and,  hence,  the  breaking 
is  suppressed  to  slightly  higher  amplitudes.  This  can  be  seen  by  examining  the 
limiting  case  of  an  “exact”  waveform  in  Fig.  3b,  where  a  ymax  —  ymin  =  -722  is 
maintained  without  breaking.  However,  a  very  small  subharmonic  disturbance 
would  cause  these  high  waves  to  break. 

Since  the  total  energy  is  constant  throughout  an  entire  numerical  simulation  as 
well  as  a  (presumably  inviscid)  experiment,  it  would  appear  to  be  an  ideal  criterion 
to  determine  breaking.  Unfortunately,  without  a  very  carefully  calibrated  and 
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instrumented  wavemaker  or  the  ability  to  measure  the  velocity  everywhere  in  the 
flow  field,  the  total  energy  cannot  be  measured.  Instead,  a  measured  steepness 
criterion  (>7max  ~  >7min)/waveiength  is  usually  used.  We  see  from  Fig.  5a  that 
this  criterion  varies  widely  in  time  for  the  two  cases  with  initial  conditions  (4) 
of  a  =  0.27  and  0.28.  The  peak-to-peak  hdght  for  the  nonbreaking  (nonspilling) 
wave  goes  higher  than  the  value  at  a  previous  time  for  a  wave  that  breaks.  Since 
(especially  the  piecewise-linear)  computations  waves  show  the  wave  breaking  at 
less  than  the  maximum  peak-to-peak  height,  the  height  for  a  nonbreaking  wave 
could  exceed  that  of  a  wave  that  is  breaking.  There  is  experimental  evidence  for 
this  as  well  (Melville  and  Rapp,  1988). 

However,  the  potential  energy,  although  not  constant  in  time,  is  much  less  vari¬ 
able  and  can  easily  be  determined  from  wave  probe  data.  Fig.  5b  demonstrates 
that  the  computed  potential  energy  for  these  same  two  initial  conditions  are  dis¬ 
tinct,  indicating  that  potential  energy  is  a  better  criterion  in  determining  whether 
a  traveling  wave  will  break.  These  computations  show  that  breaking  does  not  occur 
at  the  peak  of  the  potential  energy.  This  could  be  anticipated  since  the  increased 
fluid  velocities  near  the  crest  would  increase  the  kinetic  energy  at  the  expense  of 
the  potential  energy. 

It  should  be  noted  that  the  maximum  and  minimum  wave  heights  are  deter¬ 
mined  at  the  same  instant  of  time  at  two  different  locations  for  the  numerical 
results  of  Fig.  5a,  while  experimental  measurements  (with  one  wave  probe)  are 
usually  measured  at  one  location  at  two  different  times. 

Figs.  5a  and  5b  also  compare  spectral  and  piecewise-linear  computations.  The 


20 


piecewise-linear  computations  proceed  further  before  breaking  down,  but  as  these 
computations  sire  refined,  they  approach  the  spectral  computations  and  do  not 
proceed  as  far. 

4.4  Convergent  wave  channel 

The  growth  rate,  7,  of  equation  (7)  will  cause  a  wave  of  any  nonzero  initial  ampli¬ 
tude  to  break  eventually.  When  the  growth  rate  is  large,  the  wave  quickly  plunges 
(Fig.  6a);  when  it  is  smaller,  after  a  longo^  time  the  wave  spills  like  those  seen 
experimentally  (Fig.  6b).  The  time  required  for  the  wave  to  break,  of  course,  also 
depends  on  the  initial  conditions. 

Fig.  7  shows  the  temporal  development  of  ymax  —  yminy  the  potential  energy, 
and  the  total  energy  for  a  growth  parameter  7  =  0.2  and  two  different  initial 
conditions,  a  =  0.1  and  a  =  0.2.  The  average  growth  in  the  total  energy  for 
all  cases  is  exponential  at  the  rate  expected,  but  with  small  oscillations.  These 
oscillations  axe  not  computational  errors,  but  artifacts  of  the  growth  model.  The 
potential  energy  (and  hence  the  kinetic  energy)  also  grow  exponentiaJly  but  with 
larger  oscillations.  The  lines  cease  at  the  point  where  the  computation  fails,  which 
for  these  cases  result  in  wave  profiles  that  appear  to  be  spilling  breakers.  All 
computations  are  spectral  except  for  the  dotted  line,  which  shows  small  deviations 
of  the  peak-to-peak  measurement  near  breaking  for  a  =  0.2  when  the  piecewise- 
linear  algorithm  is  used.  In  contrast  to  the  no-growth  breaking  (Figs.  5a  and  b), 
these  numerical  simulations  show  breaking  at  or  near  the  maximum  values  of  the 
peak-to-peak  or  potential  energy  meaisures.  The  wave  grows  to  a  higher  value 
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(by  either  measure)  before  breaking  when  the  initial  amplitude  is  smaller.  This 
appears  to  be  caused  by  the  ability  of  the  smaller  wave  to  evolve  in  time  with  a 
wave  profile  that  is  more  similar  to  a  nonlinear  steadily  progressive  wave. 

Breaking  wave  criteria  for  waves  with  varying  growth  rates  can  be  obtained  from 
Fig.  8.  This  shows  that  the  breaking  criteria  (by  either  peak-to-peak  or  potential 
energy  measures)  increase  with  the  energy  input  rate,  except  for  a  “resonance” 
phenomena  around  7  =  0.1,  which  appears  to  be  an  artifact  of  the  periodic  con¬ 
straints.  The  data  also  show  dependence  on  the  initial  conditions — smaller  initial 
amplitudes  adjust  to  the  growth  “better”  and,  hence,  break  at  higher  amplitudes. 

The  dependence  of  the  breaking  height  has  been  correlated  with  growth  rate  by 
Van  Dom  and  Pazan  (1975),  who  have  also  conducted  experiments  in  a  convergent 
channel.  They  obtain  somewhat  higher  values  of  breaking  wave  height  than  those 
of  Ramberg  and  Griffin  (1987),  which  is  consistent  with  the  higher  convergence 
rate  of  their  channel,  e  =  1/10. 

4.5  Modelling  periodic  submerged  disturbances 

Duncan  (1981,  1983)  towed  hydrofoils  to  create  breaking  waves.  There  are  two 
reasons  for  using  hydrofoils:  the  first  is  to  avoid  separation  (which  we  does  not 
concern  our  potential  flow  model),  and  the  second  is  to  apply  (negative)  values  of 
lift.  We  could  model  lift  by  adding  a  periodic  vortex  array  to  the  periodic  dipole 
array,  but  we  choose  not  to  include  this  effect  here.  Rather,  we  “tow”  simple 
dipoles  at  different  speeds  and  depths  and  determine  the  strength  of  the  dipole 
(the  approximate  radius  of  the  cylinder)  that  causes  breaking. 
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We  tow  the  dipole  at  constant  speed  starting  from  rest  with  homogeneous  initial 
conditions.  To  avoid  an  impulse  at  t  =  0,  which  would  cause  high-frequency  waves, 
we  increase  the  dipole  strength  exponentially  with  a  time  constant  of  unity,  i.e, 

adp  =  ^Vdp[rdp(l.  -  (25) 

FVom  steady  linear  theory,  the  number  of  waves  that  should  appear  in  a  compu¬ 
tational  period  is  equal  to  l/t>3p  (for  orir  scaling).  We  have  chosen  Vdp  values  of  1, 
0.707,  and  0.5  to  represent  1,  2,  and  4  wave  computations. 

Figs.  9a-c  show  a  typical  computation  for  a  dipole  starting  at  the  origin  with 
a  speed  equal  to  0.5  and  unit  depth;  the  strength  Odp  represents  a  cylinder  of 
approximately  circular  cross  section  with  a  radius  =  0.39  (a  larger  radius  would 
cause  wave  breaking).  Fig.  9a  shows  the  development  of  the  free  surface  beginning 
at  rest;  initially,  a  single  peak  occurs  slightly  ahead  of  the  dipole.  A  local  maximum 
in  wave  height  is  a^jiieved  around  t  =  6,  and  then  the  wave  puts  energy  back  into 
the  dipole.  Fig.  9b  shows  computations  at  a  later  time  with  a  transition  between 
three  local  maxima  and  four  maxima  (four  is  predicted  by  linear,  steady  theory). 
This  transition  between  the  number  of  peaks  is  partly  responsible  for  the  “beating” 
rhythm,  as  seen  in  the  wave  diagnostics  in  Fig.  9c. 

Table  I  shows  the  maximum  value  of  adp  before  breaking  occurs  (up  to  /  =  200) 
for  three  combinations  of  dipole  depth  and  speed  as  computed  by  the  piecewise- 
linear  algorithm.  Often,  however,  the  wave  breaks  at  a  lower  height  as  in  the 
constant  energy  computations  in  Fig.  5b.  The  maximum  Odp  values  for  the  spectral 
algorithm  are  somewhat  lower,  although  the  wave  breaks  at  approximately  the 
same  height.  The  difference  between  the  algorithms  are  greater  for  the  dipole  case 
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because  the  computation  is  more  sensitive  to  the  numerical  parameters.  The  time 
step  had  to  be  significantly  smaller  in  the  computations  with  dipole  disturbances 
although  a  spatial  resolution  of  JV  =  64  stiU  seemed  sufficient. 


The  potential  and  total  energy  are  computed  as  before.  While  the  potential 
energy  representation  is  the  same  in  this  case,  the  kinetic  energy  should  exclude 
the  kinetic  energy  inside  the  cylinder  formed  by  the  dipole.  This  could  be  computed 
by  integrating  /  4>dtl>  in  the  clockwise  direction  about  the  exact  cylinder  surface. 
Then  conservation  of  energy  could  be  obtained  by  knowing  the  work  done  by  the 
dipole  using  Lagally’s  theorem.  This  was  not  attempted  here;  instead,  we  compute 
the  “total”  energy  from  the  free  surface  contour  only.  Mziss  is  still  conserved  to  a 
high  degree  (the  mean  height  is  10~^  to  10“*°  of  ♦he  RMS  wave  height). 


Table  I.  Wave  Breaking  Caused  by  a  Periodic  Dipole  Array 


speed 

1 

.707 

.5 

depth 

.5 

max  adp 

.14 

.10 

Umax  ymin 

.86 

.21 

Total  Energy 

.43 

mM 

.02 

Potential  Energy 

.21 

.035 

.009 

Beat  Period 

125 

130 

oo 

1 

max  Cdp 

.19 

.14 

KIH 

ymax  ymin 

.86 

.35 

Total  Energy 

.43 

.07 

.015 

Potential  Energy 

.21 

.03 

.008 

Beat  Period 

130 

160 

12 

2 

max  Odp 

KTiTil 

ymax  ymin 

mm 

Total  Energy 

.41 

.04 

Potential  Energy 

.20 

.022 

.015 

Beat  Period 

123 

200 

10 
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Figs.  lOa-b  show  breaking  waves  for  two  different  dipole  depths.  The  deeper  dipoles 
tend  to  produce  only  one  primary  peak  at  vjp  =  0.5,  while  shallow  dipoles  give 
four  as  predicted  from  linear  theory.  Surprisingly,  the  peak-to-peak  breaking  wave 
heights  are  about  the  same,  independent  of  the  size  of  the  other  nonbreaking  peaks. 

4.6  Wave  modulation  and  reflection 

Experimental  waves  break  at  significantly  less  than  the  predicted  Stokes  limiting 
value,  and  even  less  than  that  predicted  here,  because  of  three-dimensional  effects, 
viscous  effects,  wave  reflections  from  beaches,  and  the  Benjamin-Feir  instability. 
We  can  test  the  last  effect  by  taking  a  larger  periodic  domain  and  applpng  a  sub¬ 
harmonic  disturbance — a  modulated  initial  condition.  Dold  and  Peregrine  (1986) 
computed  recurrence  caused  by  modulation,  but  did  not  examine  breaking  waves. 
Fig.  11  shows  a  breaking  wave  that  has  a  10  percent  modulation  of  the  primary 
wave  for  a  wavelength  twice  that  of  the  primary  wave.  This  is  the  easiest  long 
wave  disturbance  to  model  using  a  periodic  algorithm  and  is  the  wavenumber  that 
Longuet-Hig^ns  (1978)  considers  to  be  most  unstable.  Specifically,  the  boimdary 
conditions  we  apply  are  those  of  (4)  except  that  a  is  a  “slowly”  varying  function 
of  X.  For  Fig.  11,  a  =  a'(l  -}-  6cos2x),  where  a!  =  .23  and  6  =  .10.  F.om  this  sim¬ 
ple  example,  we  see  that  a  10  percent  modulation  of  the  initial  conditions  causes 
an  approximately  20  percent  reduction  in  the  initi2d  aize  of  a  wave  that  breaks 
(o  =  .28  is  reduced  to  .23). 

We  also  expect  that  reflection  from  the  beadi  could  cause  smaller  amplitude 
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waves  to  break.  These  reflections  are  modeled  by 


1/ =  a  (sin  X  — /i  cos  x)  and  ^  =  a  (1  + /<)co8x,  (26) 

where  n  is  the  reflection  coeflScient.  Fig.  12  shows  that,  contrary  to  the  the  effect 
of  modulation,  a  10  percent  reflection  coefficient  modifles  the  initial  conditions  for 
breaking  by  approximately  5  percent  (a  =  .28  is  reduced  to  .265).  The  plots  of 
ymax—ymin  sJ^d  potential  energy  are  qualitatively  similar  to  those  cf  Figs.  5a  and  5b. 

5  Convergent  Channel  Experiments  Compared  to  Com¬ 
putations 

The  Naval  Research  Laboratory  (NRL)  experiments  discussed  in  this  section  were 
conducted  in  a  channel  30m  long  and  1.3m  wide  at  the  wavemaker  with  about 
1.0m  mean  water  depth.  The  test  procedures  and  the  set-up  of  the  exi>eriments 
are  described  by  Ramberg  et  al.  (1985)  and  Ramberg  and  Griffin  (1987).  The 
channel  was  fitted  with  a  convergent  section  with  a  rate  of  1:16.  This  generated 
steep  2uid  breaking  waves  under  reasonably  well-controlled  conditions.  Previous 
breaking  wave  experiments  in  a  wave  channel  with  a  convergent  section  heui  been 
conducted  and  reported  by  Van  Dom  and  Pazan  (1975). 

Wave  heights  were  measured  with  capacitance  wave  probes.  The  length-scale 
for  the  experiments  was  g(T/2iry,  such  that  the  nondimensional  spatial  period  of 
a  linear  wave  was  27r,  in  accordance  with  the  periodic  computations. 

The  location  of  wave  breaking  was  established  visually  during  the  experiments 
as  the  position  along  the  convergent  channel  where  the  sharp  wave  crests  were 
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first  perceived  to  be  "tripping”  into  a  spilling  or  plunging  mode  (where  we  first 
observed  an  increase  in  the  crest  fluid  velocity  over  that  of  the  traveling  waveform. 
These  locations  were  recorded  and  later  compared  to  the  positions  where  the  mea¬ 
sured  variation  of  the  average  wave  height  H  exhibited  a  transition  from  growth 
to  attenuation  along  the  channel.  In  all  of  the  cases  compared,  there  was  good 
agreement  between  the  two  estimated  locations.  Typical  examples  of  spilling  and 
plunging  breakers  are  given  in  the  photographs  in  Ramberg  and  Griffin  (1987). 

Two  hundred  equally  spaced  temporal  measurements  of  wave  height  were  taken 
over  eight  wavemaker  periods  with  two  probes  placed  3.0m  apart.  These  measure¬ 
ments  were  repeated  sequentially  over  twenty-four  spatial  locations. 

It  is  often  difficult  to  determine  the  breaking  location,  especially  when  breaking 
is  intermittent.  Fig.  13  shows  how  the  spectra  at  the  break  point  can  help  deter¬ 
mine  occurrences  of  wave  breaking.  The  nonbreaking  wave  has  one  prominent  peak 
at  Ar  =  8  since  data  were  taken  for  eight  Qrcles  of  the  wavemaker  at  each  location. 
When  breaking  occurs,  the  higher  harmonics  are  much  more  pronounced.  Also, 
the  wave  numbers  below  k  =  S  are  larger — indicating  greater  wave  modulation 
during  breaking  in  these  examples. 

At  each  location,  for  each  wavemaker  period,  a  maximum  and  a  minimum  height 
were  determined  using  quadratic  interpolation.  The  mean  and  standard  deviation 
of  the  nondimensional  peak-to-peak  height,  /f,  was  determined  for  the  eight  qrcles 
at  each  location.  The  potential  energy  (or  RMS  wave  height)  of  the  experimental 
waves  was  computed  by  integrating  the  square  of  the  height  using  Simpson's  rule  for 
six  to  eight  full  wave  periods,  as  determined  by  consecutive  crossings  of  the  average 
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datum  height  level.  After  computing  the  mean  potential  energy,  the  standard 
deviation  of  potential  energy  between  full  wave  periods  was  computed.  The  average 
RMS  value  for  all  data  (breaking  and  nonbreaking)  is  shown  in  Fig.  14  as  a  function 
of  the  nondimensional  growth  rate.  The  growth  rate  was  determined  from  (9)  even 
though  some  data  was  taken  just  into  the  nonconvergent  portion  of  the  channel  as 
described  by  Ramberg,  et  al.  (1985).  Those  waves  determined  to  be  at  the  point 
of  incipient  breaking  are  marked  with  closed  symbols.  The  general  trend  of  higher 
waves  for  higher  growth  rates  is  an  indication  of  the  growth  of  the  waves  down  the 
channel  as  the  width  decreases. 

Schultz  et  al.  (1986)  and  Ramberg  and  Griffin  (1987)  found  that  the  dimensional 
peak-to-peak  measurement  for  breaking  waves  was  equal  to  0.021  correspond¬ 
ing  to  H  =  .83  in  the  present  scaling  as  noted  in  Fig.  15a.  This  figure  shows  that 
H  increases  with  the  growth  parameter  7,  as  predicted  in  Fig.  7.  These  predictions 
are  higher  than  the  experimental  data — presumably  due  to  three-dimoisional,  vis¬ 
cous,  and  wave  reflection  and  wave  modulation  effects  in  the  experiments.  The 
least-squares  linear  fit  of  the  breaking  data  in  Fig.  15a  is  given  by 

^  =  0.57-1-7.67.  (27) 

The  standard  deviation  about  this  curve  is  ±10  percent  compared  to  ±17  percent 
for  the  entire  data  breaking  data  of  Fig.  15a.  The  corresponding  RMS  height 
(proportional  to  the  square  root  of  potential  energy)  is  shown  for  breaking  waves 
in  Fig.  15b.  A  least-squares  straight  line  fit  for  this  data  has  the  form 

RM5  =  0.22 +  1.37,  (28) 
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with  a  deviation  of  db7  percent,  significantly  better  than  the  variation  for  H.  Also, 
potential  energy  (or  RMS)  reduces  scatter  in  the  standard  deviation  of  the  mea¬ 
surement  at  one  location.  At  each  breaking  event  location,  standard  deviations  of 
H  and  the  RMS  wave  height  for  the  six  to  eight  cycles  are  calculated.  The  mean 
of  these  standard  deviations  is  3.2  percent  for  H  =  ymax  —  ymin  and  2.7  percent 
for  the  RMS  height.  The  computations  (e.g..  Figs.  6b,  9c,  12)  also  show  that 
potential  energy  (or  equivalently  RMS  wave  haght)  is  “better  behaved”  than  the 
peak-to-peak  values. 

Some  comparisons  of  the  experimental  wave  breaking  location  showed  good 
agreement  (within  30  percent)  of  the  computed  time  of  breaking  multiplied  by  the 
group  velocity,  when  a  was  chosen  &om  linear  steady  wavemaker  theory. 

6  Momentum  and  Energy  Considerations  After  Breaking 

A  wide  range  of  wave  heights  at  the  onset  of  breaking  is  shown  in  Ramberg  and 
Griffin  (1987)  for  a  number  of  recent  investigations  corresponding  to  deep  water 
waves.  The  data  from  all  of  the  experiments  employ  a  crest-to-preceding-trough 
value  for  H.  The  wave  heights  measured  by  Ochi  and  Tsai  (1983)  for  the  breaking 
of  steep  nonlinear  waves  in  a  uniform  channel  cover  the  range  of  gT^  =  200  to  800 
cm,  while  the  wave  heights  measured  in  the  NRL  experiments  cover  the  range  of 
gT^  =  550  to  1100  cm.  The  wave  heights  measured  by  Duncan  (1983)  are  in  the 
range  of  gT'^  =  100  to  400  cm  and  those  measured  by  Bonmarin  and  Ramamon- 
jiarisoa  (1985)  are  in  the  range  of  gT"^  =  350  to  650  cm.  The  latter  experiments 
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were  conducted  to  measure  the  breaking  of  waves  in  a  uniform  channel,  while  in  the 
experiments  of  Duncan  the  breaking  waves  were  generated  by  towing  a  hydrofoil 
through  still  water  at  various  submergence  depths.  Above  a  critical  submergence 
depth  of  the  hydrofoil,  wave  breaking  occurred  spontaneously  (Duncan,  1983). 

It  is  also  possible  to  compare  the  results  of  Ramberg  and  Griffin  (1987)  with 
the  measurements  of  momentum  flux  losses  by  Melville  and  Rapp  (1985).  The 
latter  reported  measxirements  of  wave  breaking  in  packets  of  waves  generated  in 
a  uniform  channel  by  the  superposition  of  Fourier  components  over  a  small  band 
centered  on  the  frequency  /<>.  The  integrated  wave  amplitude  was  assumed  by 
Melville  and  Rapp  to  be  a  measure  of  the  momentum  flux,  AS,  due  to  breaking 
and  was  determined  locally  by  taking  the  difference  between  the  incipient  breaking 
condition  and  conditions  farther  downstream  (normalized  by  the  upstream  refer¬ 
ence  value  So)-  The  rate  of  momentum  flux  loss  can  be  expressed  as  (Melville 
and  Rapp,  1985) 

,  _  A  S/Sq 

^  AA:„(x-Xfr)’ 

where  ko  =  Wo/fl'  deep  water  waves,  when  Uo  is  the  angular  frequency.  A 
comparison  between  the  experiments  of  Ramberg  and  Griffin  and  those  of  Melville 
and  Rapp  (1985)  can  be  made  directly  from  the  results  in  Table  II.  Only  the  initial 
rates  of  energy  loss  are  compared  here  because  during  the  convergent  channel 
experiments  the  energy  and  momentum  losses  rarely  reached  an  asymptotic  value, 
i.e.,  the  waves  continued  to  break  throughout  the  measurement  interval.  It  should 
be  noted  that  the  rate  of  momentiun  flux  loss  from  the  experiments  of  Ramberg 
and  Griffin  is  derived  from  estimates  of  potential  energy  based  upon  measured 


30 


wave  heights.  However,  there  is  only  a  difference  of  2t  between  the  loss  rates  of 
potential  energy  and  momentum  flux. 

The  results  of  the  two  experiments  are  in  agreement  that  the  initial  rates  of 
momentum  flux  loss  due  to  plunging  breakers  are  typically  about  twice  those  of 
spilling  breakers.  Moreover,  the  results  of  Melville  and  Rapp  ako  indicate  that 
the  total  losses  for  the  two  types  of  breakers  differ  by  about  the  same  factor.  The 
smaller  momentum  flux  losses  reported  by  Melville  and  Rapp  are  expected  since 
their  results  are  for  the  integral  momentum  flux  of  a  group  of  waves.  It  is  possible 
to  qualitatively  reconcile  the  loss  rates  magnitudes  given  in  Table  II  if  about  half 
of  the  most  energetic  waves  in  each  packet  are  actually  breaking  at  any  given  time. 


Table  II.  Momentum  Flux  Losses  in  Wave  Breaking 


Breaker  Type 

Wave  Steepness 

ak 

Momentum  Flux  Loss, 

Spilling 

0.30-^ 

0.31-0.38-^-^ 

0.023+  0.048++ 

Plunging 

0.39 

0.33-0.38 

0.045+  0.089++ 

■‘‘Melville  and  Rapp  (1985),  ■’■'•■Ramberg  and  GriflBn  (1987). 


7  Concluding  Remarks 

Computations  indicate  that  the  potential  energy  of  siurface  gravity  waves  an'e  a 
better  criterion  for  the  onset  of  breaking  for  steep  nonlinear  waves  than  the  wave 
height.  The  computed  wave  height  or  steepness  appear  to  have  more  erratic  varia¬ 
tions  in  time  than  the  potential  energy.  The  potential  energy  also  reduces  the  ap¬ 
parent  scatter  observed  in  laboratory  measurements  of  individual  breaking  events. 
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There  axe  three  experimental  indications  that  the  square  root  of  the  potential 
energy  is  better  than  the  peak-to-peak  (or  steepness)  criteria  in  predicting  the 
onset  of  breaking  events.  These  are  (1)  a  better  correlation  to  a  least-squares  fit 
with  X  with  only  two-thirds  of  the  scatter  that  the  peak-to-peak  correlation  gives, 
(2)  a  smaller  percentage  standard  deviation  around  the  mean  value  at  an  individual 
location  (breaking  or  nonbreaking),  and  (3)  a  smaller  mean  percentage  variation 
of  individual  breaking  events  from  a  breaking  criteria. 

Either  brealdng  criteria  (peak-to-peak  or  potential  energy)  have  been  shown  to 
depend  on  the  energy  input  rate,  with  the  smallest  vcdues  occurring  when  the  en¬ 
ergy  input  rate  is  small  and  spilling  breakers  are  expected.  This  dependence  further 
explains  the  scatter  of  breadcing  criteria  in  the  convergent  channel  experiments. 

Continuing  studies  show  that  the  peak-to-peak  wave  heights  and  potential  en¬ 
ergy  that  can  be  sustained  without  breaddng  axe  relatively  independent  of  the 
method  of  wave  formation.  We  find  that  the  computations  tend  to  show  a  higher 
breaking  criteria  than  the  experiments,  although  this  difference  can  be  lessened 
somewhat  by  modelling  wave  modulation  and  to  a  lesser  degree,  wave  reflection. 

The  work  at  The  University  of  Michigan  was  supported  by  Naval  Research 
Laboratory  Contract  No.  N00014-85-K-Z019,  ONR  Ocean  Engineering  Division 
Contract  N00014-87-0509,  and  the  Program  in  Ship  Hydrodynamics  at  The  Uni¬ 
versity  of  Michigan,  funded  by  the  University  Research  Initiative  of  the  OflBce 
of  Naval  Research,  Contr2M:t  No.  N000184-86-K-0684.  The  authors  acknowledge 
D.T.  Rowley’s  contribution  to  the  data  reduction  and  L.  Pall’s  editorial  assistance. 
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Figure  Captions 


1  Problem  Dommn 

a.  physical  (()  space 

b.  mapped  (C)  space 

2  Comparison  of  errors  for  energy  conservation  {N  =  32) 

-  plecewise-linear  (total  cpu  seconds  =  186) 

-  spectroJ,  high  time  resolution  (total  cpu  seconds  =  28) 

—  •  —  •  —  •  —  spectral,  low  time  resolution  (total  cpu  seconds  =  100) 

3  Steadily  progressing  wave  profiles 
3'*  J/mox  ilmin  —  0.1.  (spectral) 

b.  ymax  —  ymin  =  0.115.  (piecewisc-linear) 

4  Free  surface  profiles  (At  =  A,  N  =  32) 

a.  spilling  breaker  (a  =  0.3) 

b.  plunging  breaker  (a  =  0.544) 


5a  Peak-to-peak  wave  heights  for  a  —  0.27,  0.28 

- spectral  N  =  32, - piecewise-linear  AT  =  32, . N  =  64 

5b  Potential  energy  for  same  two  initial  conditions 

- spectral  N  =  32, - piecewise-linear  N  =  32, . N  =  64 


6a  Wave  profiles  with  rapid  exponential  growth 

7  =  .5,  increment  between  plotted  profiles  At  =  0.2 
(except  last  increment  At  =  0.1) 

6b  Wave  profiles  with  slow  exponential  growth 

7  =  .02,  increment  between  plotted  profiles  At  =  0.2 

7  Evolution  of  wave  diagnostics  with  growth  parameter  7  =  0.2 
two  initial  wave  heights,  a  =  .1  ,  .2,  spectral  computations 

ymax  ymin 

-  potential  energy 

—  •  —  •  —  •  total  energy 

.  piecewise-linear  results  for  ymax  —  ymin 

8  Breaking  criteria  for  exponential  growth  conditions 
0  —  a  =  .l  A  —  a  =  .2 

ymax  ymin 

-  potential  energy 
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9a  Wave  profiles  caused  by  a  moving  dipole  (initial  devdopment) 

Vdp  =  *5,  rdp  =  .39,  ddp=  1,  0  >  t  >  10 
9b  Wave  profiles  caused  by  a  moving  dipole  (intermediate  development) 

Vdp  -  .5,  rdp  =  .39,  ddp  =  1,  20  >  <  >  23 
9c  Wave  diagnostic  for  profiles  caused  by  a  moving  dipole 

ymax  ymin 

-  potential  energy 

—  .  —  .  —  •  —  total  energy 

.  corresponding  curves  for  a  breaking  wave  rjp  =  .42 

10a  Breaking  profile  for  shallow  dipole  (v^p  =  .5,  rjp  =  .11,  ddp  =  .5) 

10b  Breaking  profile  for  deep  dipole  (vjp  =  .5,  rjp  =  .99,  djp  =  2.0) 

11  Breaking  modulated  periodic  wave 

initial  modulation  is  10  percent  with  wavelength  twice  that  of  the  primary  wave. 
-  initial  wave  profile 

-  profile  at  breaking,  t  =97.0,  97.2,  97.4,  97.6 

12  Breaking  wave  with  reflection 

reflection  coefficient  is  10  percent  with  initial  amphtude  a  =  .265 

ymax  ymin 

-  potential  energy 

13  Fourier  coefficients  for  wave  data  (experiments  from  Ramberg  et  al.,  1985) 
-  breaking  wave 

-  nonbreaking  wave 

14  RMS  wave  height  for  all  data  (experiments  from  Ramberg  et  al.,  1985) 

(solid  symbols  denote  incipient  breaking  waves) 

15a  Peak-to-peak  height  for  breaking  waves  (experiments  from  Ramberg  et  al.,  1985) 
15b  RMS  height  results  for  breaking  waves  (experiments  from  Ramberg  et  al.,  1985) 
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Fig.  9b 
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